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p I The theoretical framework of cosmology is mainly defined by gravity, of 

^ ■ which general relativity is the current model. Recent tests of general relativity 
within the A Cold Dark Matter (CDM) model have found a concordance between 
^ , predictions and the observations of the growth rate and clustering of the cosmic 
webMJ'^. General relativity has not hitherto been tested on cosmological scales 
Q I independent of the assumptions of the ACDM model. Here we report observation 
^ ■ of the gravitational redshift of light coming from galaxies in clusters at the 99 per 
cent confidence level, based upon archival data^. The measurement agrees with 
p |, the predictions of general relativity and its modification created to explain cosmic 
Q ! acceleration without the need for dark energy (/(-R) theorj^, but is inconsistent 
^ ' with alternative models designed to avoid the presence of dark matter'^ 1^. 

^, 

According to the theory of general relativitj^, light emitted from galaxies moving in the 
gravitational potential well of galaxy clusters is expected to be redshifted proportionally to the 
difference in gravitational potential $ between the clusters and an observer, i.e., Zgj. = A$/c^, 
. where c is the velocity of light in vacuum. For typical cluster masses of ~ IO^^Mq, where 
' Mq is the Sun's mass, the gravitational redshift is estimated to bd^'^'^ cZg,- ~ 10 km s~^ 



which is around two orders magnitude smaller than the Doppler shift due to the random 
O ' motions of galaxies in clusters. The method of disentangling the kinematic Doppler effect from 
gravitational redshift relies on the fact that the former gives rise to a symmetric broadening 
of the observed velocity distribution, whereas the latter shifts its centroid. A critical factor 
in detecting such a velocity shift is the number of galaxies with spectroscopically measured 
^ I velocities and the number of galaxy clusters. Both should be sufficiently high in order to 
\ reduce the error due to the Doppler width of the velocity distribution and eliminate the 
sensitivity to irregularities in cluster structure, e.g. substructures, asphericity. 

The data are compiled from the SDS^ Data Release 7 and the associated Gaussian 
Mixture Brightest Cluster Galaxy catalogu^^ containing the positions and redshifts of galaxy 
clusters identified in the survey. The cluster sample is richness-limited with a threshold 
corresponding to a cluster mass of IO^^Mq. The mean, 5- and 95-percentile values of the 
cluster richnes^^ are 16, 8, and 86 and correspond to cluster masses of around 2 x lO^^M©, 
IO^^Mq and 10^^ Mq. The typical number of spectroscopic redshifts per cluster (within a 6 
Mpc aperture and a ±4000 km s~^ velocity range around the mean cluster velocities) varies 
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from 10 for low-richness clusters to 140 for the richest ones. 

Fig. 1 shows the histograms of galaxy velocities calculated in four bins of the projected 
cluster-centric distance centred at 0.6, 1.6, 3.3 and 5.2 Mpc. The cluster centres and redshifts 
were approximated by the coordinates and redshifts of the brightest cluster galaxies, hereafter 
BCGs. The observed velocity distributions consist of two clearly distinct parts: a quasi-flat 
distribution of galaxies not belonging to the clusters (observed due to projection effect) and 
a quasi- Gaussian component associated with galaxies gravitationally bound to the cluster^. 
The latter is expected to reveal the signature of gravitational redshift in terms of a systematic 
shift of its velocity centroid. Analysis of mock kinematic data generated from cosmological 
simulations shows that the number of redshifts and clusters is sufficient to reduce all expected 
sources of noise such as substructures, cluster asphericity, non-negligible off-set between BCGs 
and clusters centres^ (both in the position on the sky and redshift space), and to allow for 
detection of gravitational redshift at nearly 3a confidence level (see SI). 

We search for gravitational redshift by measuring the mean velocity A of the quasi- 
Gaussian component of the observed velocity distribution. We carry out a Monte Carlo 
Markov Chain analysis of the data using a two-component model for the velocity distribution 
which includes both a contribution from the cluster and non-cluster galaxies (SI). Constraints 
on the mean velocity are obtained by marginalising the likelihood function over the set of 
nuisance parameters defining the shape of both components of the velocity distribution. The 
best fitting models of the velocity distributions are shown in Fig. 1 and the resulting mea- 
surements of the mean velocity as a function of the projected cluster-centric distance R are 
presented in Fig. 2. The obtained mean velocity is negative at all radii with a clear tendency 
to decline with increasing radius. The negative values arise from the fact that the rest frames 
of the clusters are defined by the observed velocities of the central galaxies. This choice of the 
reference frame implies that the gravitational redshift manifests itself ClS db blueshiflJISl (nega- 
tive mean velocity) varying with the projected cluster-centric distance from at the cluster 
centre to — |$(0)|/c at large projected radii R. 

The detection of gravitational redshift is significant at the 99 per cent confidence level. 
The integrated signal within the 6 Mpc aperture amounts to A = —7.7 ± 3.0 km s^^ which 
is consistent with the gravitational potential depths of simulated galaxy clusters ol^ A = 
— (5 — 10) km s~^. A more quantitative comparison with theoretical predictions requires 
explicit information about the mean gravitational potential profile and the distribution of 
cluster masses in the sample. We make use of the velocity dispersion profile of the composite 
cluster to constrain both functions. Then we calculate the gravitational redshift in terms 
of the mean velocity A by convolving the individual profiles of the clusters with their mass 
distribution (SI). The resulting profile (red profile in Fig. 2; see also discussion on the effect of 
the anisotropy of galaxy orbits in SI) is fully consistent with the gravitational redshift inferred 
from the velocity distributions. The fact that the same gravitational potential underlies galaxy 
motions and gravitational redshift of photons in clusters provides observational evidence of 
the equivalence principle on the scale of galaxy clusters. 
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We confront the obtained constraints on gravitational redshift with the predictions of 
alternative theories of gravity. We consider two popular models of gravity, the tensor-vector- 
scalar (hereafter TeVeS) theorj^^l^ and the f{R) modeP, designed to alleviate the problem 
of dark matter or to recover the expansion history of the Universe, respectively. Theoretical 
profiles of gravitational redshift are calculated using the relations between the generalised 
gravitational potentials of these models and the Newtonian potential (SI). The Newtonian 
potential is inferred from the observed velocity dispersion profile of the composite cluster 
under the assumption of the most reliable anisotropic model of galaxy orbits (see SI for more 
details), and constitutes the reference basis for the calculations. For TeVeS we assume that 
the total masses of galaxy clusters make up 80 per cent of those recovered under assumption of 
the Newtonian gravity. This factor lowers the ratio of dynamical-to-baryonic mass in galaxy 
clusters to the value resulting from fitting Modified Newtonian Dynamics^- (to which TeVeS 
is a relativistic generalisation) to cluster dataP^. The resulting profile of gravitational redshift 
does not match the data, deviating from the observations at the 95 percent confidence level 
(the blue dashed line in Fig. 3). This discrepancy increases with projected radius and is mostly 
caused by a logarithmic divergence of the scalar field in the regime of small accelerations, i.e., 
g < Qq and Oq ~ 10~^° m s~^, which is responsible for a 1/r modification of the gravitational 
acceleration. This result points to a critical problem for TeVeS (or Modified Newtonian 
Dynamics) in recovering the true gravitational potential at large distances around the cluster 
centres. Considering the f{R) model, we choose the least favourable set of free parameters 
maximising the departure from Newtonian gravitational acceleratiorP'. Despite this choice, 
the resulting profile of gravitational redshift is consistent with the data (the blue solid line in 
Fig. 2). 

The obtained constraints on gravity are consistent with recent tests verifying the con- 
cordance between gravity, cosmological model and observations of the large scale structure 
of the Univers^'^. An important advantage of using gravitational redshift effect is that this 
method does not depend on cosmology (see also SI) allowing to probe gravity in a direct way. 
In particular, this implies that the discrepancy between TeVeS theory and the observation^ 
is unlikely to be a consequence of a specific choice of cosmological parameters, but indeed 
points to the inadequacy of this model to describe the Universe on very large scales. 

Our results complement a series of experiments and observations aimed at confirmation 
of the predicted gravitational redshift on different scales of the Universe. Fig. 3 shows a 
summary of the detections in terms of the relative accuracy of the measurements as a function 
of the scale of the gravitational potential well. The positions of data points vary from 20 m, 
for the first ground-based experimentP^'^, to the 1 — 10 Mpc scale for galaxy clusters. On the 
scale of the solar radius we plot the measurement of gravitational redshift for the Sun^, and 
on the 2 orders of magnitude smaller scale constraints from the observations of Sirius B white 
dwarPl 1201 and space-borne hydrogen masei'^. These results make gravitational redshift the 
only effect predicted by general relativity which has been confirmed on spatial scales spanning 
22 orders of magnitude. Studying this effect in more detail relies on the size of the redshift 
sample and therefore will be possible with the advent of the next generation redshifts surveys, 
e.g., the EUCLID satellite. 
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Figure 1 Velocity distributions of galaxies combined from 7,800 SDSS galaxy clusters. The line- 
of-sight velocity (vios) distributions are plotted in four bins of the projected cluster-centric distances 
R. They are sorted from the top to bottom according to the order of radial bins indicated in the 
upper left corner and offset vertically by an arbitrary amount for presentation purposes. Red lines 
present the histograms of the observed galaxy velocities in the cluster rest frame and black solid 
lines show the best fitting models. The model assumes a linear contribution from the galaxies 
which do not belong to the cluster and a quasi-Gaussian contribution from the cluster members 
(see SI for more details). The cluster rest frames and centres are defined by the redshifts and the 
positions of the brightest cluster galaxies. The error bars represent Poisson noise. 
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Figure 2 Constraints on gravitational redshift in galaxy clusters. The effect manifests itself as a 
blueshift A of the velocity distributions of cluster galaxies in the rest frame of their BCGs. Velocity 
shifts were estimated as the mean velocity of a quasi-Gaussian component of the observed velocity 
distributions (see Fig. 1). The error bars represent the range of A parameter containing 68 per cent 
of the marginal probability and the dispersion of the projected radii in a given bin. The blueshift 
(black points) varies with the projected radius R and its value at large radii indicates the mean 
gravitational potential depth in galaxy clusters. The red profile represents theoretical predictions of 
general relativity calculated on the basis of the mean cluster gravitational potential inferred from 
fitting the velocity dispersion profile under the assumption of the most reliable anisotropic model 
of galaxy orbits (see SI for more details). Its width shows the range of A containing 68 per cent 
of the marginal probability. The blue solid and dashed lines show the profiles corresponding to two 
modifications of standard gravity: f{R) theory^ and the tensor-vector-scalar (TeVeS) modeP'^. 
Both profiles were calculated on the basis of the corresponding modified gravitational potentials 
(see SI for more details). The prediction for f{R) represents the case which maximises the deviation 
from the gravitational acceleration in standard gravity on the scales of galaxy clusters. Assuming 
isotropic orbits in fitting the velocity dispersion profile lowers the mean gravitational depth of the 
clusters by 20 per cent. The resulting profiles of gravitational redshift for general relativity and 
f{R) theory are still consistent with the data and the discrepancy between prediction of TeVeS 
and the measurements remains nearly the same. The arrows show characteristic scales related to 
the mean radius of the virialized parts of the clusters. 



8 



1.6 - 
1.4 - 
1.2 
1 

0.8 
0.6 
0.4 



of # 



Zg in clusters/GR 



o 
□ 



Zg in clusters/f(R) 
Zg in clusters/TeVeS 

y/GR+ACDIVI^ 
Eg/GR+ACDiyi^ 







<5' 



[i 



10 



■15 



10 



■10 



10"^ 10 
scale [pc] 







10^ 



Figure 3 Tlie measured-to-predicted ratio of tlie gravitational redsliift. The figure shows the 
results of different observations or experiments as a function of the spatial scale of the gravitational 
potential well. Blue and red symbols refer to detections of gravitational redshift Zg in: ground- 
based experiment!^ (blue circle), observations of Sirius B white dwari'^^ ^° (blue triangles), space- 
based experiment!^ (blue square), observation of the Sun^^ (blue star), analysis of the cluster 
data reported in this work (red circle). All measurements are compared with the predictions of 
general relativity (solid symbols). Results obtained for galaxy clusters are also compared with the 
predictions of /(i?) theory and TeVeS model (red empty symbols). As a measure of gravitational 
redshift in galaxy clusters we used the signal integrated within the aperture of 6 Mpc. The green 
square and circle show the measurement of the rate of growth of cosmic structure^ and the probe 
of gravity Eg combining the properties of galaxy-galaxy lensing, galaxy clustering and galaxy 
velocities. Both results are compared with the prediction of general relativity with a standard 
ACDM cosmological model. All error bars represent standard deviations. The relative accuracy of 
the measurement from space-born experiments^ is beyond the resolution of the plot and amounts 
to 10"^ 
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Supplementary Information 



Measuring gravitational redshift of galaxies in clusters 

The gravitational potential depth in typical galaxy clusters, expressed in terms of the velocity 
shift, is estimated at arouncP^ —10 km s^^. This is two orders of magnitude smaller than 
the Doppler shift arising from random motions of galaxies in clusters. Bearing in mind that 
the signature of the gravitational redshift lies in the mean of the velocity distribution one 
can show that in order to reduce the error on the gravitational redshift to the level of the 
effect itself one needs at least 10^ velocities. This number should grow by a factor of 4, if 
one requires a 2a detection of the effect, and probably by another factor of 2 in order to 
account for a non-negligible number of background galaxies which do not contribute to the 
effect, but give rise to the shape of the observed velocity distribution. Needles to say, the 
only means to collect as many as 10^ — 10^ galaxy velocities in clusters is stacking redshift 
data of sufficiently large number of clusters (> 10^ clusters, assuming that current redshift 
surveys typically provide 10 redshifts per cluster). 

Combining redshift data from a number of clusters also allows to reduce the error result- 
ing from local irregularities of velocity distributions in individual clusters. Such irregularities 
arise naturally from the presence of substructures or filaments along the line of sight, devia- 
tion from spherical symmetry, residual streaming motions etc. In order to address the impact 
of these factors on the error of the gravitational redshift estimate, one needs to refer to cosmo- 
logical simulations. Such an analysis was carried out by Kim & Croftf^ who concluded that 
the minimum number of galaxy clusters required to confirm the gravitational redshift effect 
at the 2(7 confidence level is ~ 3000. With this number of clusters the gravitational redshift 
may be traced up to 6 Mpc, which is 3 — 4 times larger than the size of the virialized part 
of clusters- a natural boundary condition for all methods of the mass measurement based on 
the assumption of virial equilibrium. 

Another important source of inaccuracy in the measurement of gravitational redshift, 
which was partly taken into account by Kim & Croft^°, is the choice of a central cluster 
galaxy which ideally would be an object at rest at the bottom of the gravitational potential 
well. In our work, we approximate such galaxy by a brightest cluster galaxy, hereafter BCG. 
In general, such choice is not fully justified because the positions and velocities of BCGs 
exhibit some deviations from those defined by the cluster mass centre^^'^. For example, 
the typical dispersion of the random velocities of BCGs may reach 30 — 40 per cent of the 
total velocity dispersion in galaxy cluster^. Yet, among all cluster galaxies, BCGs are those 
whose positions and velocities coincide mostly with the location and velocities of the cluster 
centres. In order to reduce the error caused by the non-vanishing velocities of BCGs to the 
level required for detection of the gravitational redshift effect, one needs to combine the data 
from a sufficiently large number of clusters, e.g., around 2500 clusters for a 2a detectioiJ^. 
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Data 



In order to compile statistically uniform and possibly the largest sample of galaxy redshifts 
in clusters, we make use of the SDS^ Data Release 7 whose integral part is a flux-limited 
spectroscopic survey providing the redshifts of nearly million galaxies brighter than Petrosian 
r-magnitude 17.77 over the area 7400 deg^. The positions and redshifts of galaxy clusters come 
from a Gaussian Mixture Brightest Cluster Galaxy cluster catalogu^^ which is the most up- 
to-date catalogue of galaxy clusters assembled on the basis of the SDSS DR7. The catalogue 
comprises 55,000 galaxy clusters at redshifts 0.1 < z < 0.55 detected by means of searching 
for red-sequence galaxies and BCGs. It provides positions and redshifts of BCGs residing in 
galaxy clusters selected from the SDSS above a certain richness limit. For the purpose of 
our analysis, we neglect all clusters whose BCGs do not have spectroscopic redshifts. This 
reduces the cluster sample by 63 per cent without affecting the relative fractions of poor and 
rich clusters. 



We approximate the coordinates and redshift of cluster centres by the positions and 
redshifts of BCGs. Then we search for all galaxies within a 6 Mpc aperture around the cluster 
centres. The radius of this aperture is ~ 3.5 larger then the virial radius and corresponds 
to the turn-around radiu^ at which the expansion of the Universe starts to dominate over 
peculiar velocities of galaxies and the velocity cut-off separates all potential cluster galaxies 
from the galaxies of background or foreground. In order to separate potential cluster members 
from distant interlopers (galaxies of background or foreground), we select only those galaxies 
whose velocities v\os in the rest frame of a related BCG, i.e. 

V\os = C— , (1) 

-1- + zbcg 

where z and ^^bcg are the redshifts of a given galaxy and related BCG respectively, he within 
the ±4000 km s~^ range. This velocity cut-off is sufficiently wide to include all cluster mem- 
bers with no respect to the cluster mass (the minimum velocity cut-off corresponding to 
the most massive clusters at small radii is around ±3000 km s~^; however, wider velocity 
range is required for precise modelling of the interloper contribution to the observed velocity 
distribution^^). 



Our final sample comprises 7, 800 clusters with the mean redshift of 0.24 and on average 
16 galaxies with spectroscopic redshift per cluster. Those clusters with less than 5 redshifts 
were not included into the sample. As the final step we combine redshift data of all clusters 
into one. The velocity diagram of the resulting composite cluster is shown in Supplementary 
Fig. 1. 
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Analysis of the velocity distributions 



To place constraints on the mean of the observed velocity distribution of cluster galaxies in- 
duced by the gravitational redshift effect we carry out a Monte Carlo Markov Chain (MCMC) 
analysis of the velocity distribution with the likelihood function defined by 

N 

L = l[f{v,^,,,\A,ci), (2) 

1=1 

where /(fios|A,a) is a model of the velocity distribution, A is the mean of the velocity 
distribution of cluster galaxies, a is a vector of nuisance parameters describing the shape of 
the velocity distribution and is the number of redshifts. In order to account for the presence 
of the interlopers (galaxies of foreground or background observed due to projection effect), 
we use the following two-component model of the velocity distributiorP' 

f{v\os) = (1 - Pci)/b(^^ios|a) + pcifci{vios\A, a), (3) 

where /ci(fios) and /b(fios) are the velocity distributions of the cluster members and the 
interlopers, respectively (both normalised to 1), and pci is a free parameter describing the 
probability of a given galaxy to be a cluster member. The choice of the functional form 
of fh{v) depends on the operational definition of cluster membership and may vary from 
a wide Gaussian distributiorP^, if one regards all galaxies beyond the virial sphere as the 
background, to a uniform distributiorJ^, if only gravitationally unbound galaxies contribute 
to the background. For the purpose of our study a uniform background is appropriate, since all 
gravitationally bound galaxies, regardless of their positions with respect to the virial sphere, 
contribute to the expected signal of gravitational redshift. However, close inspection of the 
data reveals that a uniform model of the background must be generalised to account for a 
subtle asymmetry between the number of interlopers with negative and positive velocities. 
This asymmetry arises from the fact that a flux-limited limited survey tends to include slightly 
more galaxies which are closer and, therefore, have negative Hubble velocities in the cluster 
rest frame. We find that in order to account for this effect it is sufficient to assume that 
/b(fios) is linear in velocity fios- The slope of this velocity distribution is the second nuisance 
parameter of our model, after p^. We also note that the velocity distribution of cluster 
members cannot exhibit a similar asymmetry because it is dominated by random motions of 
galaxies which are independent of the positions with respect to the cluster centre. 

The observed velocity distributions of cluster galaxies are not Gaussian at all projected 
radii R (see Fig. 2). This deviation from Gaussianity is expected and arises mostly from 
combining data from clusters of different massed, from the fact that the radial bins are 
wider than the scale of variation of the velocity dispersion profile and from an intrinsic non- 
Gaussianity of velocity distributions of individual clusters^^. Modelling these effects is beyond 
the scope of this work and, for our purpose, it is sufficient to invoke a phenomenological model 
of /ci(^^ios) providing a satisfactory fit to the data. We find that approximating /ci(t^ios) by a 
sum of two Gaussians with the same mean velocity A satisfies this condition. This introduces 
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three additional parameters into the model given by eq. ([3]): two velocity dispersions and the 
ratio of the relative weights of both Gaussian components. Performing a K-S test, we verified 
that the final fits of our model are fully consistent with the data {p = 0.99). 

We carry out the MCMC analysis of the velocity distributions in 4 radial bins of the 
projected cluster-centric distance (see Fig. 2). The number of redshifts in subsequent bins 
varies from 15, 000 in the two innermost bins to 45, 000 for the remaining two. The choice of 
these numbers is motivated by finding a balance between bin spacing and the local number 
of cluster galaxies (proportional to pd which varies from 0.9 at the cluster centre to 0.3 at 
-R > 3 Mpc). BCGs were not included in the first bin, otherwise the estimate of the mean 
would be biased towards 0. The number of clusters contributing to the subsequent bins varies 
from ~ 1000 in the innermost bin to ~ 2000 in the outermost one. 

For the conversion between the angular and physical scales we adopted a flat ACDM 
cosmology with fim = 0.3 and the Hubble constant Hq = 70 km s~^ Mpc~^. We note, however, 
that all galaxy clusters used in this work lie at low redshifts [z ^ 0.2) and, therefore, the 
impact of using a particular cosmological model on the final results is negligible. 



Gravitational redshift profile 



Assuming spherical symmetry, the gravitational redshift profile of a single galaxy cluster (in 
terms of velocity blueshift of the velocity distribution) can be calculated using the following 
formulaP 

where R is the projected cluster-centric distance, $(r) is the gravitational potential, p(r) and 
S(i?) are the 3D and surface (2D) density profiles of galaxies. In order to estimate this effect 
for the data combined from a cluster sample, one needs to convolve this expression with the 
distribution of cluster masses in the sample. Then the resulting profile of the blueshift takes 
the following form 

/A3(i?)S(i^)(dN/dM0dM. 
^ ' /S(i?)(diV/dM,)dM, ' ^ ' 

where is the virial mass and dN/ dM^ is the mass distribution. We note that Ag depends 
implicitly on the virial mass and the shape of the gravitational potential. The virial mass 
and radius are defined in terms of the overdensity parameter 6c = ?>M^/ [Aur^pc), where pc 
is the present critical density. In our calculations we adopted 5c = 102 (see e.g. Lokas & 
HoffmaiP). 

The main unknown factor in equation ([5]) is the mass distribution. This may be esti- 
mated by means of dynamical modelling of the observed velocity dispersion profile. The left 
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panel of Supplementary Fig. 2 shows the velocity dispersion profile estimated in radial bins 
(black points) by fitting eq. with /ci(f) approximated by a single Gaussian. The profile 
is truncated at R = 1.2 Mpc which is the virial radius corresponding to the anticipated lower 
limit of all virial masses in the sample which is around lO^^M©. We note that the velocity dis- 
persion profile is flatter than typical profiles observed in single galaxy clusters. This property 
arises naturally from the fact that the cluster sample is not uniform in terms of the cluster 
mass (more massive clusters give rise to growth of velocity dispersions at large radii). Here 
we use this effect to place constraints on the mass distribution in the cluster sample. 



In analogy with equation (E]), one can show that the velocity dispersion profile for 
kinematic data combined from a set of clusters can be expressed as 

_ f JdURMR)idN/dM^)dM^ ^ 1/2 
a,os[n)-^ /S(i?)(dN/dM.)dM. J ' ^""^ 

where aios(^) refers to the velocity dispersion profile of a single cluster. We approximate 
the mass distribution by a power-law, i.e., dN/dM^ oc M~°, where a is a free parameter. 
This parameterisation is mostly motivated by the fact that it resembles the observed cluster 
counts as a function of the virial mas^. In order to account for the richness threshold of 
the cluster catalogue and cosmological decay of the mass function at high masses we impose 
cut-offs on the mass distribution at low and high masses respectively. The cut-offs are fixed 
at IO^^Mq and 2 x lO^^M© which are the limits of the mass range spanned by the clusters 
of the maxBCG catalogu^^'^'^^ the predecessor of the Gaussian Mixture Brightest Cluster 
Galaxy catalogu^ used in this work. 



In order to calculate the velocity dispersion profile of a single cluster, aios(-R) in eq. Qj, 
we make use of a model of the distribution function presented by Wojtak et alW^. The model 
is constructed under assumption of spherical symmetry, constant mass-to-light ratio (galaxies 
trace dark matter) and for a wide range of possible profiles of the orbital anisotropy. We 
approximate the dark matter density profile by the NFW formulaP^, i.e. p(r) oc l/[(r/rv)(l -|- 
Cyrv)^], where Cv is the concentration parameter. Since fitting the velocity dispersion profile 
does not allow to constrain the mass profile and the anisotropy of galaxy orbits at the same 
time (the problem known as the mass-anisotropy degeneracy), we fix all parameters related to 
the orbital anisotropy. In order not to loose generality of our analysis we consider two models 
of the orbital anisotropy: anisotropic- with the anisotropy parameter P{r) = 1 — o'g{r)/a'^{r) 
varying from in the centre to 0.4 at the virial radius, and isotropic with /3(r) = 0, where cxr 
and ae are the radial and tangential velocity dispersions. We note that these two models are 
two limiting cases of a whole family of the anisotropy profiles found both in simulation^^lIM] 
and observationJ^I^ 



We evaluate the velocity dispersion numerically as the second moment of the projected 
phase-space density at fixed projected radii R. Then we correct all resulting dispersions for 
the effect of non- vanishing random velocities of BCG^. Such correction relies on replacing all 
velocity dispersions (Jios(-R) by (crios(-R)^ + cr|cQ)^/^, where ctbcg is a typical velocity dispersion 
for BCGs. For our analysis, we assume that ctbcg equals 35 per cent of the total velocity 
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dispersion within the virial radius. 

Fitting the velocity dispersion profile to the data, we obtain constraints on the concen- 
tration parameter and the slope of the mass distribution (the right panel of Supplementary 
Fig. 2). For simplicity, we assume that the concentration parameter does not vary with the 
virial mass. The obtained concentrations are smaller by 20 — 30 per cent than those found 
for simulated dark matter halos-®. This effective flattening of the mass profile most (smaller 
concentration parameters) arises from the fact that the gravitational potential of clusters is 
measured with respect to BCGs whose positions exhibit random off-sets from the true cluster 
centreJ^. 

We make use of the constraints on and a to calculate the profile of gravitational 
redshift given by eq. (E]), where the gravitational potential takes the NFW form, i.e. 

$(r) = -(GM.K)V2^(c.)^/^i^^il±^ (7) 

and 5'(cv) = l/[ln(l + Cv) — Cv/ (1 + Cv)], and the number density of galaxies is proportional to 
that of dark matter. The resulting A profiles are shown in Supplementary Fig. 3. The blue 
(red) profile corresponds to an isotropic (anisotropic) model of galaxy orbits and the widths 
of the profiles are the la ranges obtained by marginalising over all free parameters of the 
model. Both profiles are fully consistent with the constraints on A inferred from the observed 
velocity distributions (black points). Since the systematic errors induced by a choice of the 
anisotropy model are much smaller than the random errors associated with the gravitational 
redshift measurement, in the final comparison we only consider the profile calculated for an 
anisotropic model of galaxy orbits- a more reliable model in recovering the kinematics of 
galaxy cluster^. We also note the gravitational potential resulting from the analysis with 
anisotropic galaxy orbits is the basis for the calculations of the gravitational redshift effect in 
alternative models of gravity (see the last section of SI). 

The theoretical calculations presented in this section rely on the extrapolation of the 
NFW profile beyond the virial sphere which is still robust to 2rv, but is probably less justified 
for larger radii^. In order to check the potential impact of this assumption on the estimation 
of A(i?), we considered a set of simple modifications of the density profile at large radii. 
We found that varying the asymptotic slopes of the dark matter (or galaxy number) density 
profile by ±0.5 at r > Ir^ does not change the A profiles by more than 4 per cent. We also 
found that the mass limits imposed on the mass distribution have a negligible effect on the 
final A profile. We verified that changing both mass limits by 50 per cent induces error which 
do not exceed the uncertainty due to the unknown anisotropy of galaxy orbits. 
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Test on a mock data sample 

In order to test the statistical robustness of the gravitational redshift detection, we analyse 
mock kinematic data generated from cosmological A^-body simulations of a standard ACDM 
cosmological model (for details of the simulations see Wojtak et alW^. Such test allows to 
check whether all effects related to the internal substructure of clusters and their perturbed 
velocity distributions are sufficiently reduced in the procedure of stacking kinematic data of 
a large number of galaxy clusters. 

We generate mock redshift data by drawing randomly dark matter particles from 7, 800 
cylinders of observation (the number of cylinders corresponds to the number of clusters in 
the real data sample compiled from the SDSS). The cylinders of observations are defined by 
the ±4000 km s~^ velocity range and the 6 Mpc aperture. Their viewing angles are chosen 
at random and geometrical centres are located at the centre of one of 80 cluster-mass dark 
matter halos. In order to simulate the offset between the positions of the cluster central 
galaxies and the cluster mass centres, we introduce a random shift between the positions and 
velocities of the cylinders and the corresponding halo mass centres. The maximum value of 
this shift is 0.1 Mpc for the positions— and ±35 per cent of the total velocity dispersion 
for velocities projected onto the line of sightP^. The total number of velocities in a composite 
cluster is fixed at 10^ or 10^. The former corresponds to the number of redshifts in the current 
SDSS sample and the latter represents a forecast for the future. 

In order to perform a test of detectability of the gravitational redshift effect, we consider 
two data samples. The first consists of velocities which are not corrected for the gravitational 
redshift and the second takes into account an additional velocity shift proportional to the 
local gravitational potential. In both cases particle positions remain the same. Both mock 
data samples are analysed in the same manner as the SDSS data. The obtained constraints 
on the mean velocity as a function of the projected radius are shown in Supplementary Fig. 4, 
where blue and red colour refer to the samples with 10^ and 10^ redshifts, respectively. 

The mean velocity obtained for the sample without gravitational redshift (points with 
the dashed error bars) are consistent with A = at all radii. This shows that the number of 
clusters used in our work is sufficient to reduce all local effects giving rise to the fluctuations 
of the mean velocity. Velocity shifts obtained for the second data sample (points with the 
solid error bars) clearly indicate the presence of the gravitational redshift effect. They also 
trace the true proflle of this effect averaged over the halo sample used for generating mock 
kinematic data (black solid line). We flnd that the measured velocity shifts for 10^ data 
points (blue points with the solid error bars) deviate from A = proflle at nearly the same 
signiflcance level as the results obtained for the SDSS data. We also note that an increase of 
the number of redshifts per cluster improves the errors of the measurement according to the 
rule of the inverse square root proportionality. 
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Gravitational redshift in alternative models of gravity 

The only difference in the calculation of the gravitational redshift for models of modified 
gravity relies on replacing the Newtonian gravitational potential in eq. (jlj) by the potential 
emerging from a given gravity. For the f{R) model, the gravitational potential is calculated 
using phenomenological relations obtained by Schmidlf^ who ran a series of cosmological 
simulations of this model and quantified its effect on cluster dynamics in terms of the effective 
enhancement of gravitational acceleration with respect to the standard gravity. In our work we 
consider the most critical case of his results (with |d//di?| ~ 10""^) leading to a homogeneous 
amplification of gravitational acceleration at all radii by factor of 1.33. 

Gravitational potential in TeVeS theory is a sum of the Newtonian potential $n and 
the scalar field 0. The former is calculated assuming 80 per cent of the total mass inferred 
from the velocity dispersion profile in the framework of Newtonian gravity. This factor lowers 
the total-to-barynic mass ratio to the value estimated in galaxy clusters under assumption of 
the Modified Newtonian Dynamic^. Assuming spherical symmetry, the scalar field is a 
solution of the following equatiorP 

MV4> = (fc/47r)V<l>N, (8) 

where y = {'V(f))'^kP, k << 1 and \/k/{A7Tl) ^ ao = 10"^^ m s~^. In our calculations we 
adopted fi{y) = y/y/{l + y/y), which is one of the commonly used interpolating function in 
Modified Newtonian Dynamicd^, and k = 0.01. 

In both schemes of the calculation, the reference gravitational potential of Newto- 
nian gravity is given by the constraints from fitting the velocity dispersion profile with an 
anisotropic model of galaxy orbits. We note that using the results of an isotropic model does 
not change the main conclusion of the test summarised in Fig. 3: the f{R) model is still con- 
sistent with the data, whereas TeVeS yields even more divergent profile of the gravitational 
redshift. 
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Supplementary Figure 1 Velocity diagram combined from kinematic data of 7800 galaxy clusters 
detected in the SDSS^"^ Data Release 7. Velocities vios of galaxies with respect to the brightest 
cluster galaxies are plotted as a function of the projected cluster-centric distance R. Blue lines are 
the iso-density contours equally spaced in the logarithm of galaxy density in the fios — R plane. 
The arrows show characteristic scales related to the mean virial radius estimated in dynamical 
analysis of the velocity dispersion profile. Data points represent 20 per cent of the total sample. 
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Supplementary Figure 2 Velocity dispersion profile of the composite cluster (left panel) and 
constraints on the concentration parameter and the logarithmic slope of the mass distribution 
OL (right panel) from fitting the velocity dispersion profile with an isotropic (blue) or anisotropic 
(red) model of galaxy orbits. The solid lines in the left panel show the best-fitting profiles of the 
velocity dispersion profile. The contours in the right panel are the boundaries of the \o and 2a 
confidence regions of the likelihood function. The error bars in the left panel represent the range 
containing 68 per cent of the marginal probability. 
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Supplementary Figure 3 Theoretical predictions for gravitational redshift in terms of the mean 
velocity A. The blue and red profiles were calculated on the basis of the mean cluster gravitational 
potential inferred from fitting the velocity dispersion profile with an isotropic and anisotropic model 
of galaxy orbits, respectively The widths of both profiles represent the range of A containing 68 
per cent of the probability marginalised over parameters of the gravitational potential and the 
distribution of clusters masses in the sample. The arrows show characteristic scales related to the 
mean virial radius r^. Black points show constraints on gravitational redshift from the analysis of 
the observed velocity distributions. The error bars represent the range of A parameter containing 
68 per cent of the marginal probability and the dispersion of the projected radii in a given bin. 
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Supplementary Figure 4 Constraints on the velocity shift in mock composite clusters generated 
from cosmological A^-body simulations. Blue and red points correspond to the data samples with 
10^ and 10^ redshifts, respectively. Data points with the solid (dashed) error bars represent the 
results obtained for the samples which include (do not include) the effect of gravitational redshift. 
Black solid line shows the true profile of the velocity shift A caused by gravitational redshift. The 
error bars represent the range of A parameter containing 68 per cent of the marginal probability 
and the dispersion of the projected radii in a given bin. The arrows show characteristic scales 
related to the mean virial radius r^. 
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